//========================================================================================
// Mer.cc
// ----------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/src/Mer.cc#4 $
// ------------------------------------------------------------------------------
// $Log: $
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//

#include "Mer.hh"

int Mer::mMerSize = 32; // Default mer size is 32 

//----------------------------------------------------------------------------------------
// * Mer
//----------------------------------------------------------------------------------------
//
Mer::Mer(string s){
  Sequence Seq;
  Seq.InitFromString(s);
  
  if (Seq.size() != Mer::mMerSize) {
    cout<<"Error.  Attempt to init Mer with different sized sequence."<<endl;
    cout<<"Seq.size="<<Seq.size()<<" MerSize="<<Mer::mMerSize<<endl;
    return;
  }
  //	cout<<hex;

  mLower32 = 0;
  mUpper32 = 0;
  register unsigned int loc;
  register unsigned int x;
  register int MSize = Mer::mMerSize;
  for(register int i = 0;i< MSize;i++){
    //x = Base2Bit(Seq[i]);
    x= (Seq[i] & 0x06) >> 1;
    // Does this base go to mLower32 or mUpper32?
    if (i < 16){
      loc = i<<1;               // figure out where this will go in string. 		 
      x <<= loc;              // Move those bits there. 
      mLower32 = mLower32 | x;  // Will set the desired bits to the value of x

      //		cout<<"Lower: "<<Seq[i]<<"\t"<<x<<"\t"<<mLower32<<endl;

    }else{
      loc = i-16;
      loc = loc<<1;
      x<<=loc;                 // Move those bits there.
      mUpper32 = mUpper32 | x;  // Will set the desired bits to the value of x		  

      //		cout<<"Upper: "<<Seq[i]<<"\t"<<x<<"\t"<<mUpper32<<endl;

    }  
  }
}

//----------------------------------------------------------------------------------------
// * Mer
//----------------------------------------------------------------------------------------
// Constructor to create a 2-bit encoding of up to 32-mer.   Smaller mers, if supported, 
// will have their upper bits set to all zero.  32 mers are spread over two 32 bit words.  
// Bases 0-15 in Lower32, with base 0 at bit 0, base 1 at bit 2, etc.  Bases 16-31 in 
// Upper32 with base 16 at bit 0, base 17 at bit 2, etc.
//
Mer::Mer(Sequence &Seq){
  if (Seq.size() != Mer::mMerSize) {
    cout<<"Error.  Attempt to init Mer with different sized sequence."<<endl;
    cout<<"Seq.size="<<Seq.size()<<" MerSize="<<Mer::mMerSize<<endl;
    return;
  }

  mLower32 = 0;
  mUpper32 = 0;
  register unsigned int loc;
  register unsigned int x;
  register int MSize = Mer::mMerSize;
  for(register int i = 0;i< MSize;i++){
    //x = Base2Bit(Seq[i]);
    x= (Seq[i] & 0x06) >> 1;
    // Does this base go to mLower32 or mUpper32?
    if (i < 16){
      loc = i<<1;               // figure out where this will go in string. 		 
      x <<= loc;              // Move those bits there. 
      mLower32 = mLower32 | x;  // Will set the desired bits to the value of x
    }else{
      loc = i-16;
      loc = loc<<1;
      x<<=loc;                 // Move those bits there.
      mUpper32 = mUpper32 | x;  // Will set the desired bits to the value of x		  
    }
  
  }
}

//----------------------------------------------------------------------------------------
// * Mer
//----------------------------------------------------------------------------------------
// Constructor to create a 2-bit encoding of a specified subsequence of a sequence.   
//
Mer::Mer(Sequence3Bit &Seq,int merStart,bool &bContainsNX){

  bContainsNX = false;
  if ((Seq.mBaseSize-merStart)< Mer::mMerSize){
    cout<<"Error.  Attempt to init Mer with different sized sequence."<<endl;
    cout<<"Seq.size="<<Seq.mBaseSize<<" MerStart="<<merStart<<" MerSize="<<Mer::mMerSize<<endl;
    return;
  }
  
  mLower32 = 0;
  mUpper32 = 0;
  register unsigned int loc;
  register unsigned int x;
  register int MSize = Mer::mMerSize;
  for(register int i = 0 ;i< MSize;i++){
    // See if it contains a N or X
    // KJD BUG  Shouldn't this check for lower case also???
    if (Seq[i+merStart] & 0x04){
      bContainsNX = true;
      return;
    }

    // 0x03 masks out all but bits 0 and 1, which in 3 bit encoding is where base is. 
    x= (Seq[i+merStart] & 0x03);
    // Does this base go to mLower32 or mUpper32?
    if (i < 16){
      loc = i<<1;               // figure out where this will go in string. 		 
      x <<= loc;              // Move those bits there. 
      mLower32 = mLower32 | x;  // Will set the desired bits to the value of x
    }else{
      loc = i-16;
      loc = loc<<1;
      x<<=loc;                 // Move those bits there.
      mUpper32 = mUpper32 | x;  // Will set the desired bits to the value of x		  
    }
  } 

}

//----------------------------------------------------------------------------------------
// * Mer
//----------------------------------------------------------------------------------------
// Constructor to create a 2-bit encoding of a specified subsequence of a sequence.   
//
Mer::Mer(Sequence3Bit &Seq,int merStart){
  if ((Seq.mBaseSize-merStart)< Mer::mMerSize){
    cout<<"Error.  Attempt to init Mer with different sized sequence."<<endl;
    cout<<"Seq.size="<<Seq.mBaseSize<<"merStart="<<merStart<<" MerSize="<<Mer::mMerSize<<endl;
    return;
  }
  
  mLower32 = 0;
  mUpper32 = 0;
  register unsigned int loc;
  register unsigned int x;
  register int MSize = Mer::mMerSize;
  for(register int i = 0 ;i< MSize;i++){

    // Converting 3Bit to 2Bit just want to drop bit 2.  
    x= Seq[i+merStart] & 0x03;

    // Does this base go to mLower32 or mUpper32?
    if (i < 16){
      loc = i<<1;               // figure out where this will go in string. 		 
      x <<= loc;              // Move those bits there. 
      mLower32 = mLower32 | x;  // Will set the desired bits to the value of x
    }else{
      loc = i-16;
      loc = loc<<1;
      x<<=loc;                 // Move those bits there.
      mUpper32 = mUpper32 | x;  // Will set the desired bits to the value of x		  
    }
  } 
}

//----------------------------------------------------------------------------------------
// * Mer
//----------------------------------------------------------------------------------------
// Constructor to create a 2-bit encoding of a specified subsequence of a sequence.   
//
Mer::Mer(Sequence &Seq,int merStart){
  int diff = (Seq.size()-merStart);
  register int merSize = Mer::mMerSize; // To avoid future possible de-reference. 

  if (diff < merSize){		
    cout<<"Error.  Attempt to init Mer with different sized sequence."<<endl;
    int diff = (Seq.size()-merStart);
    cout<<"diff="<<diff<<endl;
    cout<<"Seq.size="<<Seq.size()<<"\tmerStart="<<merStart<<" MerSize="<<Mer::mMerSize<<endl;
    return;
  }
  
  mLower32 = 0;
  mUpper32 = 0;
  register unsigned int loc;
  register unsigned int x;
  for(register int i = 0 ;i< merSize;i++){
    //x = Base2Bit(Seq[i+merStart]);

    x= (Seq[i+merStart] & 0x06) >> 1;
    // Does this base go to mLower32 or mUpper32?
    if (i < 16){
      loc = i<<1;               // figure out where this will go in string. 		 
      x <<= loc;              // Move those bits there. 
      mLower32 = mLower32 | x;  // Will set the desired bits to the value of x
    }else{
      loc = i-16;
      loc = loc<<1;
      x<<=loc;                 // Move those bits there.
      mUpper32 = mUpper32 | x;  // Will set the desired bits to the value of x		  
    }
  } 
}


//----------------------------------------------------------------------------------------
// * Mer
//----------------------------------------------------------------------------------------
// Constructor to create a 2-bit encoding of a specified subsequence of a sequence.   
//
Mer::Mer(Sequence &Seq,int merStart,bool &bContainsNX){

	bContainsNX = false;

  int diff = (Seq.size()-merStart);
  register int merSize = Mer::mMerSize; // To avoid future possible de-reference. 

  if (diff < merSize){		
    cout<<"Error.  Attempt to init Mer with different sized sequence."<<endl;
    int diff = (Seq.size()-merStart);
    cout<<"diff="<<diff<<endl;
    cout<<"Seq.size="<<Seq.size()<<"\tmerStart="<<merStart<<" MerSize="<<Mer::mMerSize<<endl;
    return;
  }
  
  mLower32 = 0;
  mUpper32 = 0;
  register unsigned int loc;
  register unsigned int x;
  for(register int i = 0 ;i< merSize;i++){
    //x = Base2Bit(Seq[i+merStart]);

		// See if it contains a N or X
		char testchr = Seq[i+merStart];
    if ((toupper(testchr)=='N') || (toupper(testchr) == 'X')){
      bContainsNX = true;
      return;
    }

    x= (Seq[i+merStart] & 0x06) >> 1;
    // Does this base go to mLower32 or mUpper32?
    if (i < 16){
      loc = i<<1;               // figure out where this will go in string. 		 
      x <<= loc;              // Move those bits there. 
      mLower32 = mLower32 | x;  // Will set the desired bits to the value of x
    }else{
      loc = i-16;
      loc = loc<<1;
      x<<=loc;                 // Move those bits there.
      mUpper32 = mUpper32 | x;  // Will set the desired bits to the value of x		  
    }
  } 
}








//----------------------------------------------------------------------------------------
// * ToSequence
//----------------------------------------------------------------------------------------
// Converts mer back into a sequence.
//
void 
Mer::ToSequence(Sequence &S){
  S.resize(Mer::mMerSize);
  unsigned int lowerTemp = mLower32;
  unsigned int upperTemp = mUpper32;

  register int MSize = Mer::mMerSize;
  for(int i = 0;i< MSize;i++){
    if (i < 16){
      unsigned int x = lowerTemp & 0x00000003;
      S[i] = Bit2Base(x);
      lowerTemp = lowerTemp >> 2;
    }else{
      unsigned int x = upperTemp & 0x00000003;
      S[i] = Bit2Base(x);
      upperTemp = upperTemp >> 2;
    }
  }
}

//----------------------------------------------------------------------------------------
// * Bit2Base
//----------------------------------------------------------------------------------------
// Converts a 2-bit encoding for a single base int the ASCII character for that base. 
inline char 
Mer::Bit2Base(unsigned int x){
  switch(x){
    case 0: return('a');break;
    case 1: return('c');break;
    case 2: return('t');break;
    case 3: return('g');break;
    default: return('?');break;
  };
}
//----------------------------------------------------------------------------------------
// * Base2Bit
//----------------------------------------------------------------------------------------
// Converts the ASCII character for a base into the 2-bit encoding.  Assumes that b is one 
// of A,C,T,G.   Behavior undefined otherwise.  
//
inline unsigned int
Mer::Base2Bit(char b){
  return((b & 0x06) >> 1);
}

//----------------------------------------------------------------------------------------
// * ComplementBits
//----------------------------------------------------------------------------------------
// Computes the complement on the 2-bit representation. 
//
inline unsigned int 
Mer::ComplementBits(unsigned int bits)
{
  switch(bits){
    case 0: return(2);
    case 1: return(3);
    case 2: return(0);
    case 3: return(1);	 
  }
}


//----------------------------------------------------------------------------------------
// * ReverseComplement
//----------------------------------------------------------------------------------------
//
void
Mer::ReverseComplement(Mer &fwdMer)
{
  mUpper32 = fwdMer.mUpper32;
  mLower32 = fwdMer.mLower32;
  ReverseComplement();
}


//----------------------------------------------------------------------------------------
// * ReverseComplement
//----------------------------------------------------------------------------------------
// Reverse complement assuming bits 1/2 of ASCII as representation of a base.  It is 
// slower computing RC with this encoding than it would be encoding {A,C,G,T} = {0,1,2,3}, 
// where complement is just logical negation, but the latter encoding is slower on Mer creation
// which is the dominant time.   
//
void
Mer::ReverseComplement()
{
  register unsigned int tUpper = mUpper32;
  register unsigned int tLower = mLower32;
  register unsigned int newUpper=0,newLower =0;
	
  register int MSize = Mer::mMerSize;

  // There are two cases, all the bases fit in one word, or they take 
  // some or all of two words.  
  if (MSize <= 16){

    // Copy (with rev-comp) 
    register unsigned int shift = (2*MSize)-2;
    for(register int i = 0; i < MSize;i++){
      register unsigned int bits = tLower & 0x00000003;  // Grab two lowest bits. 
      bits = ComplementBits(bits);  // complement them. 
      bits <<=shift;     // Shift them to the right position...
      newLower = newLower | bits; // Add them to the new Lower
      tLower>>=2;
      shift-=2;   
    }
  }else {
    // Note:  overlap is in bases.  		
    register int overlap = MSize - 16;
		
    // Copy (with rev-comp) overlapping bits from lower into upper...
		
    // Note: shift is in bits, not bases. 
    register unsigned int shift = (2*overlap) - 2;
    for(register int i = 0; i < overlap;i++){
      register unsigned int bits = tLower & 0x00000003;  // Grab two lowest bits. 
      bits = ComplementBits(bits);  // complement them. 
      bits <<=shift;     // Shift them to the right position...
			
      newUpper = newUpper | bits; // Add them to the new Upper

      tLower>>=2;
      shift-=2;   
    }
		
    // Copy (with rev-comp) non-overlapping bits from lower into lower...
    register unsigned int coverlap = 16-overlap;
    shift = 30;
    for(register int i = 0; i < coverlap;i++){
      register unsigned int bits = tLower & 0x00000003;  // Grab two lowest bits. 
      bits = ComplementBits(bits);  // complement them. 
      bits <<=shift;     // Shift them to the right position...
      newLower = newLower | bits; // Add them to the new Lower

      tLower>>=2;
      shift-=2;   
    }
		
    // Copy (with rev-comp) overlapping bits from upper into lower....
    shift = (2*overlap)-2;
    for(register int i = 0; i < overlap;i++){
      register unsigned int bits = tUpper & 0x00000003;  // Grab two lowest bits. 
      bits = ComplementBits(bits);  // complement them. 
      bits <<=shift;     // Shift them to the right position...
      newLower = newLower | bits; // Add them to the new Lower

      tUpper>>=2;
      shift-=2;   
    }
  }

  mUpper32 = newUpper;
  mLower32 = newLower;
}


//---------------------------------
// * InSubset
//---------------------------------
//
bool
Mer::InSubset(unsigned int Subset)
{
  // grab 4 bits from the lower32 bits, and shift them to the 
  // left three places (i.e. multiply by 8). 
  unsigned int lowerBit = (mLower32 & 0x0000000F);
  lowerBit = lowerBit << 3;

  // grab 3 bits from the upper 32 bits and add them to the 
  unsigned int mySet = (mUpper32 & 0x00000007) +lowerBit;

  // Unfortunately, left as it is, Subset=0 means all zero bits, 
  // which translates to all A's, which skews the batch.  So 
  // need to randomize the bits a bit. Does this work, or does it
  // just change which batch contains the all a's???
  // mySet = mySet ^ 0x2D;

  if (mySet == Subset) return(true);
  else return(false);
}
